*clean_population
use "Yourpath\co-est2022-alldata.dta" 
drop npopchg2020 npopchg2021 npopchg2022 births2020 births2021 births2022 deaths2020 deaths2021 deaths2022 naturalchg2020 naturalchg2021 naturalchg2022 internationalmig2020 internationalmig2021 internationalmig2022 domesticmig2020 domesticmig2021 domesticmig2022 netmig2020 netmig2021 netmig2022 residual2020 residual2021 residual2022 gqestimatesbase2020 gqestimates2020 gqestimates2021 gqestimates2022 rbirth2021 rbirth2022 rdeath2021 rdeath2022 rnaturalchg2021 rnaturalchg2022 rinternationalmig2021 rinternationalmig2022 rdomesticmig2021 rdomesticmig2022 rnetmig2021 rnetmig2022

*rename
rename stname mst
rename ctyname mcty

*check type
describe mst mcty, short

*check and delete duplicates

duplicates report mst mcty
duplicates list mst mcty
duplicates drop mst mcty, force
duplicates report mst mcty

*save file
save "Yourpath\cleaned_population.dta", replace



use "Yourpath\merged data - weekly disease Google mobility Delphi.dta"

*create mobility
egen mobility = rowmean(retail_and_recreation_percent_ch grocery_and_pharmacy_percent_cha parks_percent_change_from_baseli transit_stations_percent_change_ workplaces_percent_change_from_b)

drop _merge_cases_mobility state val_pct_wearing_mask_5d se_pct_wearing_mask_5d sample_size_pct_wearing_mask_5d represented_pct_wearing_mask_5d val_pct_wearing_mask_5d_alt se_pct_wearing_mask_5d_alt sample_size_pct_wearing_mask_5d_ represented_pct_wearing_mask_5d_ val_pct_worried_become_ill se_pct_worried_become_ill sample_size_pct_worried_become_i represented_pct_worried_become_i val_pct_anxious_5d se_pct_anxious_5d sample_size_pct_anxious_5d represented_pct_anxious_5d val_pct_anxious_5d_alt se_pct_anxious_5d_alt sample_size_pct_anxious_5d_alt represented_pct_anxious_5d_alt val_pct_depressed_5d se_pct_depressed_5d sample_size_pct_depressed_5d represented_pct_depressed_5d val_pct_depressed_5d_alt se_pct_depressed_5d_alt sample_size_pct_depressed_5d_alt represented_pct_depressed_5d_alt val_pct_felt_isolated_5d se_pct_felt_isolated_5d sample_size_pct_felt_isolated_5d represented_pct_felt_isolated_5d val_pct_felt_isolated_5d_alt se_pct_felt_isolated_5d_alt v539 v540 val_pct_vaccinated_or_accept se_pct_vaccinated_or_accept sample_size_pct_vaccinated_or_ac represented_pct_vaccinated_or_ac val_pct_wearing_mask_7d se_pct_wearing_mask_7d sample_size_pct_wearing_mask_7d represented_pct_wearing_mask_7d val_pct_wearing_mask_7d_alt se_pct_wearing_mask_7d_alt sample_size_pct_wearing_mask_7d_ represented_pct_wearing_mask_7d_ val_pct_anxious_7d se_pct_anxious_7d sample_size_pct_anxious_7d represented_pct_anxious_7d val_pct_anxious_7d_alt se_pct_anxious_7d_alt sample_size_pct_anxious_7d_alt represented_pct_anxious_7d_alt val_pct_depressed_7d se_pct_depressed_7d sample_size_pct_depressed_7d represented_pct_depressed_7d val_pct_depressed_7d_alt se_pct_depressed_7d_alt sample_size_pct_depressed_7d_alt represented_pct_depressed_7d_alt val_pct_felt_isolated_7d se_pct_felt_isolated_7d sample_size_pct_felt_isolated_7d represented_pct_felt_isolated_7d val_pct_felt_isolated_7d_alt se_pct_felt_isolated_7d_alt v1407 v1408 val_pct_accept_vaccine_no_appoin se_pct_accept_vaccine_no_appoint sample_size_pct_accept_vaccine_n represented_pct_accept_vaccine_n val_pct_vaccinated_appointment_o se_pct_vaccinated_appointment_or sample_size_pct_vaccinated_appoi represented_pct_vaccinated_appoi val_pct_worried_catch_covid se_pct_worried_catch_covid sample_size_pct_worried_catch_co represented_pct_worried_catch_co start_date end_date _merge_cases_mobility_delphi

*rename
rename statename mst
rename county mcty

*check type
describe mst mcty, short

*check type
describe wy, short

*refine wrong data
replace epi_week = "29-Feb" if epi_week == "Feb-29"

*re-categorization
drop weeklycases weeklydeaths
 
* Check if epi_year exists
capture confirm variable epi_year

* If epi_year exists, use the replace command to update it
if _rc == 0 {
    * Extract the last two digits of the year and convert it to numeric
    gen year_suffix = real(substr(from, length(from)-1, 2))
    
    * Replace the epi_year column
    replace epi_year = 2000 + year_suffix
}
* If epi_year does not exist, create it using the gen command
else {
    * Extract the last two digits of the year and convert it to numeric
    gen year_suffix = real(substr(from, length(from)-1, 2))
    
    * Generate the epi_year column
    gen epi_year = 2000 + year_suffix
}
drop year_suffix
 
* Split the from column into day, month, and year
split from, parse("-") gen(temp)

* Convert the generated temporary variables into numeric (if necessary)
gen day = real(temp1)
gen month_str = temp2
gen year = real(temp3) + 2000

* Since Stata may not directly convert month abbreviations to numeric, we will convert manually
gen month = .
replace month = 01 if month_str == "Jan"
replace month = 02 if month_str == "Feb"
replace month = 03 if month_str == "Mar"
replace month = 04 if month_str == "Apr"
replace month = 05 if month_str == "May"
replace month = 06 if month_str == "Jun"
replace month = 07 if month_str == "Jul"
replace month = 08 if month_str == "Aug"
replace month = 09 if month_str == "Sep"
replace month = 10 if month_str == "Oct"
replace month = 11 if month_str == "Nov"
replace month = 12 if month_str == "Dec"

* Use the mdy() function to generate Stata's internal date format
gen date = mdy(month, day, year)

* Set the format for the from_date column to NNDDCCYY
format date %tdNNDDCCYY

* Check the result of the conversion
list date in 1/10

* Regenerate epi_week

rename epi_week epi_week_old
* Initialize the epi_week variable
gen epi_week = .

* Set the epidemiological week for the year 2020
local start_date_2020 = mdy(1, 4, 2020)
replace epi_week = 202001 + floor((date - `start_date_2020') / 7) if year(date) == 2020 & date >= `start_date_2020'

* Set the epidemiological week for the year 2021, and so on
local start_date_2021 = mdy(1, 3, 2021)
replace epi_week = 202101 + floor((date - `start_date_2021') / 7) if year(date) == 2021 & date >= `start_date_2021'

* Note: Based on the provided rules, you will need to manually set the epi_week for the last week of each year
replace epi_week = 202053 if date >= mdy(12, 27, 2020) & date <= mdy(1, 2, 2021)
replace epi_week = 202152 if date >= mdy(12, 26, 2021) & date <= mdy(1, 1, 2022)

* For the year 2022
local start_date_2022 = mdy(1, 2, 2022)
replace epi_week = 202201 + floor((date - `start_date_2022') / 7) if year(date) == 2022 & date >= `start_date_2022'

* Special handling for the last week of 2022
replace epi_week = 202252 if date >= mdy(12, 25, 2022) & date <= mdy(12, 31, 2022)

* Set the epidemiological week for the year 2023
local start_date_2023 = mdy(1, 1, 2023)
replace epi_week = 202301 + floor((date - `start_date_2023') / 7) if year(date) == 2023 & date >= `start_date_2023'

describe epi_year epi_week, short

sort mst mcty epi_week
duplicates report mst mcty epi_week

* description of mobility
summarize mobility, detail
gen epi_year_wy = string(epi_week)

save "Yourpath\cleaned_mobility.dta", replace


*case and death
use "Yourpath\county level case and death.dta"

*clean
*rename
rename statename mst
rename county mcty 

*
describe epi_year epi_week, short
gen epi_year_wy = string(epi_year) + string(epi_week, "%02.0f")

sort mst mcty epi_year_wy
duplicates report mst mcty epi_year_wy
duplicates list mst mcty epi_year_wy
duplicates tag mst mcty epi_year_wy, generate(dup_tag)
duplicates drop mst mcty epi_year_wy, force
duplicates report mst mcty epi_year_wy

save "Yourpath\cleaned_casedeath.dta", replace

clear
use "Yourpath\cleaned_mobility.dta"

merge m:1 mst mcty using "Yourpath\cleaned_population.dta", force

drop _merge
merge 1:1 mst mcty epi_year_wy using "Yourpath\cleaned_casedeath.dta", force

*add population to CT
* Manually set population estimates for 2020 based on the county and year, ensuring mst is "Connecticut"
replace popestimate2020 = 957050 if mcty == "Fairfield County" & epi_year == 2020 & mst == "Connecticut"
replace popestimate2020 = 898682 if mcty == "Hartford County" & epi_year == 2020 & mst == "Connecticut"
replace popestimate2020 = 184938 if mcty == "Litchfield County" & epi_year == 2020 & mst == "Connecticut"
replace popestimate2020 = 164063 if mcty == "Middlesex County" & epi_year == 2020 & mst == "Connecticut"
replace popestimate2020 = 864094 if mcty == "New Haven County" & epi_year == 2020 & mst == "Connecticut"
replace popestimate2020 = 268450 if mcty == "New London County" & epi_year == 2020 & mst == "Connecticut"
replace popestimate2020 = 149767 if mcty == "Tolland County" & epi_year == 2020 & mst == "Connecticut"
replace popestimate2020 = 116404 if mcty == "Windham County" & epi_year == 2020 & mst == "Connecticut"

* Manually set population estimates for 2021 based on the county and year, ensuring mst is "Connecticut"
replace popestimate2021 = 959768 if mcty == "Fairfield County" & epi_year == 2021 & mst == "Connecticut"
replace popestimate2021 = 896854 if mcty == "Hartford County" & epi_year == 2021 & mst == "Connecticut"
replace popestimate2021 = 185000 if mcty == "Litchfield County" & epi_year == 2021 & mst == "Connecticut"
replace popestimate2021 = 164759 if mcty == "Middlesex County" & epi_year == 2021 & mst == "Connecticut"
replace popestimate2021 = 863700 if mcty == "New Haven County" & epi_year == 2021 & mst == "Connecticut"
replace popestimate2021 = 268805 if mcty == "New London County" & epi_year == 2021 & mst == "Connecticut"
replace popestimate2021 = 150293 if mcty == "Tolland County" & epi_year == 2021 & mst == "Connecticut"
replace popestimate2021 = 116418 if mcty == "Windham County" & epi_year == 2021 & mst == "Connecticut"

* Manually set population estimates for 2022 based on the county and year, ensuring mst is "Connecticut"
replace popestimate2022 = 796253 if mcty == "Fairfield County" & epi_year == 2022 & mst == "Connecticut"
replace popestimate2022 = 878018 if mcty == "Hartford County" & epi_year == 2022 & mst == "Connecticut"
replace popestimate2022 = 132322 if mcty == "Litchfield County" & epi_year == 2022 & mst == "Connecticut"
replace popestimate2022 = 158669 if mcty == "Middlesex County" & epi_year == 2022 & mst == "Connecticut"
replace popestimate2022 = 841208 if mcty == "New Haven County" & epi_year == 2022 & mst == "Connecticut"
replace popestimate2022 = 242651 if mcty == "New London County" & epi_year == 2022 & mst == "Connecticut"
replace popestimate2022 = 129797 if mcty == "Tolland County" & epi_year == 2022 & mst == "Connecticut"
replace popestimate2022 = 118792 if mcty == "Windham County" & epi_year == 2022 & mst == "Connecticut"

* Initialize the phase variable
gen phase = .

* Set phase=2 for the period (October 1, 2020 to July 31, 2021)
replace phase = 2 if date >= mdy(10, 1, 2020) & date <= mdy(7, 31, 2021)

* Set phase=1 for the period before October 1, 2020
replace phase = 1 if date < mdy(10, 1, 2020)

* Set phase=3 for the period after July 31, 2021
replace phase = 3 if date > mdy(7, 31, 2021)

* Create a new column wcp based on the first four digits of epi_year_wy
gen wcp = cond(substr(epi_year_wy, 1, 4) == "2020", weeklycases / (popestimate2020/100000), cond(substr(epi_year_wy, 1, 4) == "2021", weeklycases / (popestimate2021/100000), cond(substr(epi_year_wy, 1, 4) == "2022", weeklycases / (popestimate2022/100000), .)))

* Create a new column wdp based on the first four digits of epi_year_wy
gen wdp = cond(substr(epi_year_wy, 1, 4) == "2020", weeklydeaths / (popestimate2020/100000), cond(substr(epi_year_wy, 1, 4) == "2021", weeklydeaths / (popestimate2021/100000), cond(substr(epi_year_wy, 1, 4) == "2022", weeklydeaths / (popestimate2022/100000), .)))

* Create log-transformed column for wcp
gen lgwcp = log(wcp+1)

* Create log-transformed column for wdp
gen lgwdp = log(wdp+1)

* Plot histograms
histogram wcp
histogram lgwcp, title("Histogram of lgwcp")

histogram wdp
histogram lgwdp, title("Histogram of lgwdp")

* Plot histograms for phase 2 data
histogram wcp if phase==2, title("Histogram of wcp for phase=2")
histogram wdp if phase==2, title("Histogram of wdp for phase=2")
histogram lgwcp if phase==2, title("Histogram of lgwcp for phase=2")
histogram lgwdp if phase==2, title("Histogram of lgwdp for phase=2")

* Count non-missing observations for mobility and lgwcp
count if !missing(mobility) & !missing(lgwcp)

* Create a helper variable to mark cases where all specified variables are non-missing
gen all_nonmissingm = !missing(mobility) & !missing(lgwcp) & !missing(lgwdp)
gen all_nonmissingp = !missing(lgwcp) & !missing(lgwdp)

* Use the summarize command to view descriptive statistics for cases where all_nonmissing is 1
summarize mobility lgwcp lgwdp if all_nonmissingm
summarize mobility lgwcp lgwdp if all_nonmissingp

* Extract the month from the date
rename month old_month
gen month = month(date)

* Initialize seasonal variables
gen spring = 0
gen summer = 0
gen fall = 0

* Set seasonal variables based on the month
replace spring = 1 if inrange(month, 3, 5)
replace summer = 1 if inrange(month, 6, 8)
replace fall = 1 if inrange(month, 9, 11)


*generate cyclic

sort fips_code date
by fips_code: gen time =_n
gen cyclic=mod(time,4)
order time cyclic, after(epi_week)

egen epi_year_wy_id = group(epi_year_wy)
sort fips_code epi_year_wy_id
duplicates report fips_code epi_year_wy_id
duplicates list fips_code epi_year_wy_id
duplicates tag fips_code epi_year_wy_id, generate(dup_tag1)
duplicates drop fips_code epi_year_wy_id, force
duplicates report fips_code epi_year_wy_id

summarize mobility lgwcp if all_nonmissingm


xtset fips_code epi_year_wy_id
save "Yourpath\062024_combineddata.dta", replace



**mobility-death
clear
*state level
use "Yourpath\062024_combineddata.dta"
xtset fips_code epi_year_wy_id
xtunitroot ips mobility, demean
xtunitroot ips lgwdp, demean

 *Create a group identifier for each unique mst value
egen group_id = group(mst)

 *Determine how many unique groups (states) there are
su group_id, meanonly
local num_unique_states = r(max)

 *Divide the states into groups of 7, calculating the total number of groups
local group_size = 7
local num_groups = ceil(`num_unique_states' / `group_size')

 *Save data for each group
forvalues i = 1/`num_groups' {
     *Reload the original dataset and recreate group_id to avoid "group_id not found" error
    if `i' > 1 {
        clear
        use "Yourpath\062024_combineddata.dta", clear
        egen group_id = group(mst)
    }

     *Calculate the group number range for each group
    local start_group_id = (`i' - 1) * `group_size' + 1
    local end_group_id = min(`start_group_id' + `group_size' - 1, `num_unique_states')

     *Keep data for the current group
    keep if group_id >= `start_group_id' & group_id <= `end_group_id'

     *Save data for the current group to a new file
    save "Yourpath\state_group_`i'.dta", replace
}

 *Clear the workspace
clear

 *Define a temporary file handle
tempname handle
postfile `handle' str100 state mmsc_aic using "mmsc_results1.dta", replace

* Process data file for i=1
local i = 1
local filename "Yourpath\\state_group_`i'.dta"
clear
use "`filename'", clear
xtset fips_code epi_year_wy_id



* Get unique states
levelsof mst, local(states)
local num_states: word count `states'
di "Number of unique states: " `num_states'

* Loop over each state
foreach s in `states' {
    local modified_s = subinstr("`s'", " ", "", .)
    local all_estimates ""

    * Loop over different values of a and b
    forvalues a=1/4 {
        forvalues b=1/8 {
            quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if mst == "`s'" & phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
            local current_est "model_`modified_s'_`a'_`b'"
            estimates store `current_est'
            local all_estimates "`all_estimates' `current_est'"
        }
    }

    * Run estat mmsc and capture MMSC-AIC values
    foreach est in `all_estimates' {
        quietly estat mmsc `est'
        if _rc == 0 {
            matrix results = r(S)
            local nrows = rowsof(results)
            local ncols = colsof(results)
            di "Processing estimates for " "`est'" ": rows = " `nrows' " cols = " `ncols'
            if `ncols' >= 5 {
                * Extract MMSC-AIC values from the second row
                forval j = 2/`nrows' {
                    local mmsc_aic_value = results[`j',5]
                    post `handle' ("`modified_s'") (`mmsc_aic_value')
                }
            }
            else {
                di "Warning: Expected column not found in results matrix for `est'. Matrix cols = " `ncols'
            }
        }
        else {
            di "Failed to execute estat mmsc for `est'"
        }
    }
}
postclose `handle'

 *Define a second temporary file handle
tempname handle
postfile `handle' str100 state mmsc_aic using "mmsc_results2.dta", replace

* Process data file for i=2
local i = 2
local filename "Yourpath\\state_group_`i'.dta"
clear
use "`filename'", clear
xtset fips_code epi_year_wy_id

* Replace the value in the mst column
replace mst = "Columbia" if mst == "District of Columbia"



* Get unique states
levelsof mst, local(states)
local num_states: word count `states'
di "Number of unique states: " `num_states'

* Loop over each state
foreach s in `states' {
    local modified_s = subinstr("`s'", " ", "", .)
    local all_estimates ""

    * Loop over different values of a and b
    forvalues a=1/4 {
        forvalues b=1/8 {
            quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if mst == "`s'" & phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
            local current_est "model_`modified_s'_`a'_`b'"
            estimates store `current_est'
            local all_estimates "`all_estimates' `current_est'"
        }
    }

    * Run estat mmsc and capture MMSC-AIC values
    foreach est in `all_estimates' {
        quietly estat mmsc `est'
        if _rc == 0 {
            matrix results = r(S)
            local nrows = rowsof(results)
            local ncols = colsof(results)
            di "Processing estimates for " "`est'" ": rows = " `nrows' " cols = " `ncols'
            if `ncols' >= 5 {
                * Extract MMSC-AIC values from the second row
                forval j = 2/`nrows' {
                    local mmsc_aic_value = results[`j',5]
                    post `handle' ("`modified_s'") (`mmsc_aic_value')
                }
            }
            else {
                di "Warning: Expected column not found in results matrix for `est'. Matrix cols = " `ncols'
            }
        }
        else {
            di "Failed to execute estat mmsc for `est'"
        }
    }
}
postclose `handle'

tempname handle
postfile `handle' str100 state mmsc_aic using "mmsc_results3.dta", replace

* Process the data file for i=3
local i = 3
local filename "Yourpath\\state_group_`i'.dta"
clear
use "`filename'", clear
xtset fips_code epi_year_wy_id



* Get unique states
levelsof mst, local(states)
local num_states: word count `states'
di "Number of unique states: " `num_states'

* Loop over each state
foreach s in `states' {
    local modified_s = subinstr("`s'", " ", "", .)
    local all_estimates ""

    * Loop over different values of a and b
    forvalues a=1/4 {
        forvalues b=1/8 {
            quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if mst == "`s'" & phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
            local current_est "model_`modified_s'_`a'_`b'"
            estimates store `current_est'
            local all_estimates "`all_estimates' `current_est'"
        }
    }

    * Run estat mmsc and capture MMSC-AIC values
    foreach est in `all_estimates' {
        quietly estat mmsc `est'
        if _rc == 0 {
            matrix results = r(S)
            local nrows = rowsof(results)
            local ncols = colsof(results)
            di "Processing estimates for " "`est'" ": rows = " `nrows' " cols = " `ncols'
            if `ncols' >= 5 {
                * Extract MMSC-AIC values from the second row
                forval j = 2/`nrows' {
                    local mmsc_aic_value = results[`j',5]
                    post `handle' ("`modified_s'") (`mmsc_aic_value')
                }
            }
            else {
                di "Warning: Expected column not found in results matrix for `est'. Matrix cols = " `ncols'
            }
        }
        else {
            di "Failed to execute estat mmsc for `est'"
        }
    }
}
postclose `handle'

* Repeat the same steps for i=4, i=5

tempname handle
postfile `handle' str100 state mmsc_aic using "mmsc_results4.dta", replace

* Process the data file for i=4
local i = 4
local filename "Yourpath\\state_group_`i'.dta"
clear
use "`filename'", clear
xtset fips_code epi_year_wy_id


* Get unique states
levelsof mst, local(states)
local num_states: word count `states'
di "Number of unique states: " `num_states'

* Loop over each state
foreach s in `states' {
    local modified_s = subinstr("`s'", " ", "", .)
    local all_estimates ""

    * Loop over different values of a and b
    forvalues a=1/4 {
        forvalues b=1/8 {
            quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if mst == "`s'" & phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
            local current_est "model_`modified_s'_`a'_`b'"
            estimates store `current_est'
            local all_estimates "`all_estimates' `current_est'"
        }
    }

    * Run estat mmsc and capture MMSC-AIC values
    foreach est in `all_estimates' {
        quietly estat mmsc `est'
        if _rc == 0 {
            matrix results = r(S)
            local nrows = rowsof(results)
            local ncols = colsof(results)
            di "Processing estimates for " "`est'" ": rows = " `nrows' " cols = " `ncols'
            if `ncols' >= 5 {
                * Extract MMSC-AIC values from the second row
                forval j = 2/`nrows' {
                    local mmsc_aic_value = results[`j',5]
                    post `handle' ("`modified_s'") (`mmsc_aic_value')
                }
            }
            else {
                di "Warning: Expected column not found in results matrix for `est'. Matrix cols = " `ncols'
            }
        }
        else {
            di "Failed to execute estat mmsc for `est'"
        }
    }
}
postclose `handle'

* Repeat the same steps for i=5

tempname handle
postfile `handle' str100 state mmsc_aic using "mmsc_results5.dta", replace

* Process the data file for i=5
local i = 5
local filename "Yourpath\\state_group_`i'.dta"
clear
use "`filename'", clear
xtset fips_code epi_year_wy_id



* Get unique states
levelsof mst, local(states)
local num_states: word count `states'
di "Number of unique states: " `num_states'

* Loop over each state
foreach s in `states' {
    local modified_s = subinstr("`s'", " ", "", .)
    local all_estimates ""

    * Loop over different values of a and b
    forvalues a=1/4 {
        forvalues b=1/8 {
            quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if mst == "`s'" & phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
            local current_est "model_`modified_s'_`a'_`b'"
            estimates store `current_est'
            local all_estimates "`all_estimates' `current_est'"
        }
    }

    * Run estat mmsc and capture MMSC-AIC values
    foreach est in `all_estimates' {
        quietly estat mmsc `est'
        if _rc == 0 {
            matrix results = r(S)
            local nrows = rowsof(results)
            local ncols = colsof(results)
            di "Processing estimates for " "`est'" ": rows = " `nrows' " cols = " `ncols'
            if `ncols' >= 5 {
                * Extract MMSC-AIC values from the second row
                forval j = 2/`nrows' {
                    local mmsc_aic_value = results[`j',5]
                    post `handle' ("`modified_s'") (`mmsc_aic_value')
                }
            }
            else {
                di "Warning: Expected column not found in results matrix for `est'. Matrix cols = " `ncols'
            }
        }
        else {
            di "Failed to execute estat mmsc for `est'"
        }
    }
}
postclose `handle'


tempname handle
postfile `handle' str100 state mmsc_aic using "mmsc_results6.dta", replace

* Process the data file for i=6
local i = 6
local filename "Yourpath\\state_group_`i'.dta"
clear
use "`filename'", clear
xtset fips_code epi_year_wy_id


drop if mst == "Puerto Rico"

* Get unique states
levelsof mst, local(states)
local num_states: word count `states'
di "Number of unique states: " `num_states'

* Loop over each state
foreach s in `states' {
    local modified_s = subinstr("`s'", " ", "", .)
    count if mst == "`s'" & phase == 2 & !missing(mobility) & !missing(lgwdp)
    di "`s' number of rows meeting conditions: " _N
    if _N == 0 {
        di "Warning: No observations meet the conditions for model estimation. Skipping `s'"
        continue
    }
    local all_estimates ""
    
    * Loop over different values of a and b
    forvalues a=1/4 {
        forvalues b=1/8 {
            quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if mst == "`s'" & phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
            local current_est "model_`modified_s'_`a'_`b'"
            estimates store `current_est'
            local all_estimates "`all_estimates' `current_est'"
        }
    }

    * Run estat mmsc and capture MMSC-AIC values
    foreach est in `all_estimates' {
        quietly estat mmsc `est'
        if _rc == 0 {
            matrix results = r(S)
            local nrows = rowsof(results)
            local ncols = colsof(results)
            di "Processing estimates for " "`est'" ": rows = " `nrows' " cols = " `ncols'
            if `ncols' >= 5 {
                * Extract MMSC-AIC values from the second row
                forval j = 2/`nrows' {
                    local mmsc_aic_value = results[`j',5]
                    post `handle' ("`modified_s'") (`mmsc_aic_value')
                }
            } 
            else {
                di "Warning: Expected column not found in results matrix for `est'. Matrix cols = " `ncols'
            }
        } 
        else {
            di "Failed to execute estat mmsc for `est'"
        }
    }
}
postclose `handle'

* Repeat the same steps for i=7

tempname handle
postfile `handle' str100 state mmsc_aic using "mmsc_results7.dta", replace

* Process the data file for i=7
local i = 7
local filename "Yourpath\\state_group_`i'.dta"
clear
use "`filename'", clear
xtset fips_code epi_year_wy_id

/* Drop rows that do not meet the conditions
drop if phase != 2 | mobility == . | lgwdp == .
di "Remaining observations after drop: " _N
*/

* Get unique states
levelsof mst, local(states)
local num_states: word count `states'
di "Number of unique states: " `num_states'

* Loop over each state
foreach s in `states' {
    local modified_s = subinstr("`s'", " ", "", .)
    local all_estimates ""
    
    * Loop over different values of a and b
    forvalues a=1/4 {
        forvalues b=1/8 {
            quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if mst == "`s'" & phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
            local current_est "model_`modified_s'_`a'_`b'"
            estimates store `current_est'
            local all_estimates "`all_estimates' `current_est'"
        }
    }

    * Run estat mmsc and capture MMSC-AIC values
    foreach est in `all_estimates' {
        quietly estat mmsc `est'
        if _rc == 0 {
            matrix results = r(S)
            local nrows = rowsof(results)
            local ncols = colsof(results)
            di "Processing estimates for " "`est'" ": rows = " `nrows' " cols = " `ncols'
            if `ncols' >= 5 {
                * Extract MMSC-AIC values from the second row
                forval j = 2/`nrows' {
                    local mmsc_aic_value = results[`j',5]
                    post `handle' ("`modified_s'") (`mmsc_aic_value')
                }
            } 
            else {
                di "Warning: Expected column not found in results matrix for `est'. Matrix cols = " `ncols'
            }
        } 
        else {
            di "Failed to execute estat mmsc for `est'"
        }
    }
}
postclose `handle'

* Repeat the same steps for i=8

tempname handle
postfile `handle' str100 state mmsc_aic using "mmsc_results8.dta", replace

* Process the data file for i=8
local i = 8
local filename "Yourpath\\state_group_`i'.dta"
clear
use "`filename'", clear
xtset fips_code epi_year_wy_id

/* Drop rows that do not meet the conditions
drop if phase != 2 | mobility == . | lgwdp == .
di "Remaining observations after drop: " _N
*/

* Get unique states
levelsof mst, local(states)
local num_states: word count `states'
di "Number of unique states: " `num_states'

* Loop over each state
foreach s in `states' {
    local modified_s = subinstr("`s'", " ", "", .)
    local all_estimates ""
    
    * Loop over different values of a and b
    forvalues a=1/4 {
        forvalues b=1/8 {
            quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if mst == "`s'" & phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
            local current_est "model_`modified_s'_`a'_`b'"
            estimates store `current_est'
            local all_estimates "`all_estimates' `current_est'"
        }
    }

    * Run estat mmsc and capture MMSC-AIC values
    foreach est in `all_estimates' {
        quietly estat mmsc `est'
        if _rc == 0 {
            matrix results = r(S)
            local nrows = rowsof(results)
            local ncols = colsof(results)
            di "Processing estimates for " "`est'" ": rows = " `nrows' " cols = " `ncols'
            if `ncols' >= 5 {
                * Extract MMSC-AIC values from the second row
                forval j = 2/`nrows' {
                    local mmsc_aic_value = results[`j',5]
                    post `handle' ("`modified_s'") (`mmsc_aic_value')
                }
            } 
            else {
                di "Warning: Expected column not found in results matrix for `est'. Matrix cols = " `ncols'
            }
        } 
        else {
            di "Failed to execute estat mmsc for `est'"
        }
    }
}
postclose `handle'


* Load dataset
use "Yourpath\\mmsc_results1.dta", clear

* Ensure the data is sorted by state and another identifier (such as time or ID)
sort state 
* Generate a sequence number within each state
by state: gen seq = _n

* Calculate mode and day
gen m = ceil(seq/8)
gen d = seq - (m - 1)*8

* Convert mode and day to string
gen str_m = string(m) + "m"
gen str_d = string(d) + "d"

* Merge mode and day strings to create the num variable
capture confirm variable num

if _rc == 0 {
    replace num = str_m + str_d
}
else { 
    gen num = str_m + str_d
}

* Drop intermediate variables
drop seq m d str_m str_d

save "Yourpath\\mmsc_results1_order.dta", replace


* Repeat the same steps for mmsc_results2.dta
use "Yourpath\\mmsc_results2.dta", clear
sort state 
by state: gen seq = _n
gen m = ceil(seq/8)
gen d = seq - (m - 1)*8
gen str_m = string(m) + "m"
gen str_d = string(d) + "d"
capture confirm variable num

if _rc == 0 {
    replace num = str_m + str_d
}
else { 
    gen num = str_m + str_d
}
drop seq m d str_m str_d
save "Yourpath\\mmsc_results2_order.dta", replace


* Repeat the same steps for mmsc_results3.dta
use "Yourpath\\mmsc_results3.dta", clear
sort state 
by state: gen seq = _n
gen m = ceil(seq/8)
gen d = seq - (m - 1)*8
gen str_m = string(m) + "m"
gen str_d = string(d) + "d"
capture confirm variable num

if _rc == 0 {
    replace num = str_m + str_d
}
else { 
    gen num = str_m + str_d
}
drop seq m d str_m str_d
save "Yourpath\\mmsc_results3_order.dta", replace


* Repeat the same steps for mmsc_results4.dta
use "Yourpath\\mmsc_results4.dta", clear
sort state 
by state: gen seq = _n
gen m = ceil(seq/8)
gen d = seq - (m - 1)*8
gen str_m = string(m) + "m"
gen str_d = string(d) + "d"
capture confirm variable num

if _rc == 0 {
    replace num = str_m + str_d
}
else { 
    gen num = str_m + str_d
}
drop seq m d str_m str_d
save "Yourpath\\mmsc_results4_order.dta", replace


* Repeat the same steps for mmsc_results5.dta
use "Yourpath\\mmsc_results5.dta", clear
sort state 
by state: gen seq = _n
gen m = ceil(seq/8)
gen d = seq - (m - 1)*8
gen str_m = string(m) + "m"
gen str_d = string(d) + "d"
capture confirm variable num

if _rc == 0 {
    replace num = str_m + str_d
}
else { 
    gen num = str_m + str_d
}
drop seq m d str_m str_d
save "Yourpath\\mmsc_results5_order.dta", replace


* Repeat the same steps for mmsc_results6.dta
use "Yourpath\\mmsc_results6.dta", clear
sort state 
by state: gen seq = _n
gen m = ceil(seq/8)
gen d = seq - (m - 1)*8
gen str_m = string(m) + "m"
gen str_d = string(d) + "d"
capture confirm variable num

if _rc == 0 {
    replace num = str_m + str_d
}
else { 
    gen num = str_m + str_d
}
drop seq m d str_m str_d
save "Yourpath\\mmsc_results6_order.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results7.dta", clear

* Ensure the data is sorted by state and another identifier (such as time or ID)
sort state 
* Generate a sequence number within each state
by state: gen seq = _n

* Calculate mode and day
gen m = ceil(seq/8)
gen d = seq - (m - 1)*8

* Convert mode and day to string
gen str_m = string(m) + "m"
gen str_d = string(d) + "d"

* Merge mode and day strings to create the num variable
capture confirm variable num

if _rc == 0 {
    replace num = str_m + str_d
}
else { 
    gen num = str_m + str_d
}

* Drop intermediate variables
drop seq m d str_m str_d

save "Yourpath\\mmsc_results7_order.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results8.dta", clear

* Ensure the data is sorted by state and another identifier (such as time or ID)
sort state 
* Generate a sequence number within each state
by state: gen seq = _n

* Calculate mode and day
gen m = ceil(seq/8)
gen d = seq - (m - 1)*8

* Convert mode and day to string
gen str_m = string(m) + "m"
gen str_d = string(d) + "d"

* Merge mode and day strings to create the num variable
capture confirm variable num

if _rc == 0 {
    replace num = str_m + str_d
}
else { 
    gen num = str_m + str_d
}

* Drop intermediate variables
drop seq m d str_m str_d

save "Yourpath\\mmsc_results8_order.dta", replace

* Check the results
list state num


* Load the dataset
use "Yourpath\\mmsc_results1_order.dta", clear

* Ensure the data is sorted by state and mmsc_aic to find the smallest mmsc_aic entry within each state
sort state mmsc_aic

* Generate row number within each state
by state: gen rownum = _n

* Filter out the smallest mmsc_aic row in each state
gen min1 = rownum <= 1

* Keep the smallest mmsc_aic row and retain only state, mmsc_aic, and num columns
keep if min1
keep state mmsc_aic num

save "Yourpath\\mmsc_results1_aic_1.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results2_order.dta", clear

* Ensure the data is sorted by state and mmsc_aic to find the smallest mmsc_aic entry within each state
sort state mmsc_aic

* Generate row number within each state
by state: gen rownum = _n

* Filter out the smallest mmsc_aic row in each state
gen min1 = rownum <= 1

* Keep the smallest mmsc_aic row and retain only state, mmsc_aic, and num columns
keep if min1
keep state mmsc_aic num

save "Yourpath\\mmsc_results2_aic_1.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results3_order.dta", clear

* Ensure the data is sorted by state and mmsc_aic to find the smallest mmsc_aic entry within each state
sort state mmsc_aic

* Generate row number within each state
by state: gen rownum = _n

* Filter out the smallest mmsc_aic row in each state
gen min1 = rownum <= 1

* Keep the smallest mmsc_aic row and retain only state, mmsc_aic, and num columns
keep if min1
keep state mmsc_aic num

save "Yourpath\\mmsc_results3_aic_1.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results4_order.dta", clear

* Ensure the data is sorted by state and mmsc_aic to find the smallest mmsc_aic entry within each state
sort state mmsc_aic

* Generate row number within each state
by state: gen rownum = _n

* Filter out the smallest mmsc_aic row in each state
gen min1 = rownum <= 1

* Keep the smallest mmsc_aic row and retain only state, mmsc_aic, and num columns
keep if min1
keep state mmsc_aic num

save "Yourpath\\mmsc_results4_aic_1.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results5_order.dta", clear

* Ensure the data is sorted by state and mmsc_aic to find the smallest mmsc_aic entry within each state
sort state mmsc_aic

* Generate row number within each state
by state: gen rownum = _n

* Filter out the smallest mmsc_aic row in each state
gen min1 = rownum <= 1

* Keep the smallest mmsc_aic row and retain only state, mmsc_aic, and num columns
keep if min1
keep state mmsc_aic num

save "Yourpath\\mmsc_results5_aic_1.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results6_order.dta", clear

* Ensure the data is sorted by state and mmsc_aic to find the smallest mmsc_aic entry within each state
sort state mmsc_aic

* Generate row number within each state
by state: gen rownum = _n

* Filter out the smallest mmsc_aic row in each state
gen min1 = rownum <= 1

* Keep the smallest mmsc_aic row and retain only state, mmsc_aic, and num columns
keep if min1
keep state mmsc_aic num

save "Yourpath\\mmsc_results6_aic_1.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results7_order.dta", clear

* Ensure the data is sorted by state and mmsc_aic to find the smallest mmsc_aic entry within each state
sort state mmsc_aic

* Generate row number within each state
by state: gen rownum = _n

* Filter out the smallest mmsc_aic row in each state
gen min1 = rownum <= 1

* Keep the smallest mmsc_aic row and retain only state, mmsc_aic, and num columns
keep if min1
keep state mmsc_aic num

save "Yourpath\\mmsc_results7_aic_1.dta", replace


* Load the dataset
use "Yourpath\\mmsc_results8_order.dta", clear

* Ensure the data is sorted by state and mmsc_aic to find the smallest mmsc_aic entry within each state
sort state mmsc_aic

* Generate row number within each state
by state: gen rownum = _n

* Filter out the smallest mmsc_aic row in each state
gen min1 = rownum <= 1

* Keep the smallest mmsc_aic row and retain only state, mmsc_aic, and num columns
keep if min1
keep state mmsc_aic num

save "Yourpath\\mmsc_results8_aic_1.dta", replace


* View the results


use "Yourpath\\mmsc_results1_aic_1.dta", clear
list state mmsc_aic num 
use "Yourpath\\mmsc_results2_aic_1.dta", clear
list state mmsc_aic num 
use "Yourpath\\mmsc_results3_aic_1.dta", clear
list state mmsc_aic num 
use "Yourpath\\mmsc_results4_aic_1.dta", clear
list state mmsc_aic num 
use "Yourpath\\mmsc_results5_aic_1.dta", clear
list state mmsc_aic num 
use "Yourpath\\mmsc_results6_aic_1.dta", clear
list state mmsc_aic num 
use "Yourpath\\mmsc_results7_aic_1.dta", clear
list state mmsc_aic num 
use "Yourpath\\mmsc_results8_aic_1.dta", clear
list state mmsc_aic num





use "Yourpath\062024_combineddata.dta" 
        capture confirm variable lgwdp
        if _rc == 0 {
            egen std_lgwdp = std(lgwdp)
        }
        capture confirm variable mobility
        if _rc == 0 {
            egen std_mobility = std(mobility)
        }
		xtset fips_code epi_year_wy_id
        * 运行 xtdpdgmm 命令

xtunitroot fisher lgwdp, pperron lags(1)

xtunitroot fisher mobility, pperron lags(1)




xtdpdgmm L(0/3).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "Alabama" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/2).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Alaska" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/2).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "Arizona" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Arkansas" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "California" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Colorado" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/2).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Connecticut" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/1).std_mobility L(1/2).std_lgwdp spring fall summer i.cyclic if mst == "Delaware" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Florida" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Georgia" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/1).std_mobility L(1/4).std_lgwdp spring fall summer i.cyclic if mst == "Hawaii" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Idaho" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Illinois" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Indiana" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "Iowa" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Kansas" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/1).std_mobility L(1/4).std_lgwdp spring fall summer i.cyclic if mst == "Kentucky" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Louisiana" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "Maine" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/2).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Maryland" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/2).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Massachusetts" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "Michigan" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Minnesota" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/1).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Mississippi" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/1).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "Missouri" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/3).std_lgwdp spring fall summer i.cyclic if mst == "Montana" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Nebraska" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/3).std_lgwdp spring fall summer i.cyclic if mst == "Nevada" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "New Hampshire" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/2).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "New Jersey" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "New Mexico" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "New York" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "North Carolina" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "North Dakota" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Ohio" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Oklahoma" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Oregon" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Pennsylvania" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Rhode Island" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/4).std_lgwdp spring fall summer i.cyclic if mst == "South Carolina" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/1).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "South Dakota" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Tennessee" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Texas" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Utah" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/2).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Vermont" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "Virginia" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/4).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Washington" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/3).std_lgwdp spring fall summer i.cyclic if mst == "West Virginia" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/3).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "Wisconsin" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

xtdpdgmm L(0/1).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Wyoming" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons


use "Yourpath\062024_combineddata.dta" 
* Create standardized variables, only if they don't already exist
capture confirm variable std_lgwdp
if _rc {
    egen std_lgwdp = std(lgwdp)
}

capture confirm variable std_mobility
if _rc {
    egen std_mobility = std(mobility)
}

* Set the panel data structure
xtset fips_code epi_year_wy_id

* Manage the selection into Excel manually, but only copy data and rearrange the format

* Import the Excel file with the a and b values
import excel "Yourpath\061324_countylevel_select a b.xlsx", firstrow clear

* Save the a and b values for the states
tempfile abvalues

* Convert the mst variable to a string
encode mst, generate(mst_string)

save "Yourpath\abvalues", replace

* Load the main dataset
use "Yourpath\061124_combineddata.dta", clear
drop _merge

* Merge the a and b values into the main dataset
merge m:1 mst using "Yourpath\abvalues"

* Check the merge result
tabulate _merge

* Drop rows with missing values
drop if missing(mst)

* Trim whitespace from all values in the mst column
gen mst_clean = trim(mst)

* Replace the original mst column
drop mst
rename mst_clean mst

* Generate a local macro list of state names
levelsof mst, local(states)

* Remove spaces from state names
local states_clean
foreach state of local states {
    local states_clean `states_clean' "`trim(state)'"
}
local states "`states_clean'"

* Create a new Excel file
export excel using "Yourpath\models.csv", replace

* Loop through each state
foreach st in `states' {
    * Extract the a and b values for the state
    local a = a[`st']
    local b = b[`st']
    
    * Check if there are observations that meet the conditions
    count if mst == "`st'" & phase == 2
    if r(N) == 0 {
        display "No data available for `st'"
    }
    else {
        * Execute the model
        xtdpdgmm L(0/`a').std_mobility L(1/`b').std_lgwdp spring fall summer i.cyclic / *
                 if mst == "`st'" & phase == 2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) ///
                 iv(spring fall summer i.cyclic, d) m(d) c nocons
        
        * Store the results in Stata's temporary file
        estimates store model_`st'
        
    }
}
















* state level
* mobility_death
use "Yourpath\062024_combineddata.dta" 
* Set the panel data structure
xtset fips_code epi_year_wy_id

local all_estimates ""

* Iterate through all possible values of a, b, and c
forvalues a=1/4 {
    forvalues b=1/8 {
                * Run regression and suppress output
                quietly xtdpdgmm L(0/`a').mobility L(1/`b').lgwdp spring fall summer i.cyclic if phase==2, gmm(L.mobility L.lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
                * Store the results
                estimates store x`a'y`b'z`c'lags
                * Add the results to the string
                local all_estimates "`all_estimates' x`a'y`b'z`c'lags"
            }
        }
* Run estat mmsc command
estat mmsc `all_estimates'

egen std_lgwdp = std(lgwdp)
egen std_mobility = std(mobility)

* CountryLevel
xtdpdgmm L(0/1).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

* Use nlcom to calculate the nonlinear combination
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp]) / (1 - (_b[L1.std_mobility]))

* CountryLevel
xtdpdgmm L(0/1).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons

* Use nlcom to calculate the nonlinear combination
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp]) / (1 - (_b[L1.std_mobility]))


 *Alabama
xtdpdgmm L(0/3).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "Alabama" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Alaska
xtdpdgmm L(0/2).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Alaska" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility]))

 *Arizona
xtdpdgmm L(0/2).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "Arizona" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility]))

 *Arkansas
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Arkansas" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *California
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "California" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Colorado
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Colorado" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Connecticut
xtdpdgmm L(0/2).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Connecticut" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility]))

 *Delaware
xtdpdgmm L(0/1).std_mobility L(1/2).std_lgwdp spring fall summer i.cyclic if mst == "Delaware" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp]) / (1 - (_b[L1.std_mobility]))

 *Florida
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Florida" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Georgia
xtdpdgmm L(0/4).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Georgia" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Hawaii
xtdpdgmm L(0/1).std_mobility L(1/4).std_lgwdp spring fall summer i.cyclic if mst == "Hawaii" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp]) / (1 - (_b[L1.std_mobility]))

 *Idaho
xtdpdgmm L(0/4).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Idaho" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Illinois
xtdpdgmm L(0/3).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Illinois" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Indiana
xtdpdgmm L(0/3).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Indiana" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Iowa
xtdpdgmm L(0/4).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "Iowa" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Kansas
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Kansas" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Kentucky
xtdpdgmm L(0/1).std_mobility L(1/4).std_lgwdp spring fall summer i.cyclic if mst == "Kentucky" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp]) / (1 - (_b[L1.std_mobility]))

 *Louisiana
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Louisiana" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Maine
xtdpdgmm L(0/4).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "Maine" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Maryland
xtdpdgmm L(0/2).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Maryland" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility]))

 *Massachusetts
xtdpdgmm L(0/2).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Massachusetts" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility]))

 *Michigan
xtdpdgmm L(0/3).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "Michigan" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Minnesota
xtdpdgmm L(0/4).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Minnesota" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Mississippi
xtdpdgmm L(0/1).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Mississippi" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility]))

 *Missouri
xtdpdgmm L(0/1).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "Missouri" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp]) / (1 - (_b[L1.std_mobility]))

 *Montana
xtdpdgmm L(0/3).std_mobility L(1/3).std_lgwdp spring fall summer i.cyclic if mst == "Montana" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Nebraska
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Nebraska" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Nevada
xtdpdgmm L(0/3).std_mobility L(1/3).std_lgwdp spring fall summer i.cyclic if mst == "Nevada" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *New Hampshire
xtdpdgmm L(0/3).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "New Hampshire" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *New Jersey
xtdpdgmm L(0/2).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "New Jersey" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility]))

 *New Mexico
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "New Mexico" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *New York
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "New York" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *North Carolina
xtdpdgmm L(0/4).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "North Carolina" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *North Dakota
xtdpdgmm L(0/3).std_mobility L(1/7).std_lgwdp spring fall summer i.cyclic if mst == "North Dakota" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Ohio
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Ohio" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Oklahoma
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Oklahoma" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Oregon
xtdpdgmm L(0/3).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Oregon" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Pennsylvania
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Pennsylvania" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Rhode Island
xtdpdgmm L(0/3).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Rhode Island" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *South Carolina
xtdpdgmm L(0/4).std_mobility L(1/4).std_lgwdp spring fall summer i.cyclic if mst == "South Carolina" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *South Dakota
xtdpdgmm L(0/1).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "South Dakota" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility]))

 *Tennessee
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Tennessee" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Texas
xtdpdgmm L(0/4).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Texas" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Utah
xtdpdgmm L(0/3).std_mobility L(1/1).std_lgwdp spring fall summer i.cyclic if mst == "Utah" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Vermont
xtdpdgmm L(0/2).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Vermont" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility]))

 *Virginia
xtdpdgmm L(0/4).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "Virginia" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *Washington
xtdpdgmm L(0/4).std_mobility L(1/5).std_lgwdp spring fall summer i.cyclic if mst == "Washington" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility] + _b[L4.std_mobility]))

 *West Virginia
xtdpdgmm L(0/3).std_mobility L(1/3).std_lgwdp spring fall summer i.cyclic if mst == "West Virginia" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Wisconsin
xtdpdgmm L(0/3).std_mobility L(1/6).std_lgwdp spring fall summer i.cyclic if mst == "Wisconsin" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp]) / (1 - (_b[L1.std_mobility] + _b[L2.std_mobility] + _b[L3.std_mobility]))

 *Wyoming
xtdpdgmm L(0/1).std_mobility L(1/8).std_lgwdp spring fall summer i.cyclic if mst == "Wyoming" & phase==2, gmm(L.std_mobility L.std_lgwdp, l(1 20)) iv(spring fall summer i.cyclic, d) m(d) c nocons
nlcom (_b[L1.std_lgwdp] + _b[L2.std_lgwdp] + _b[L3.std_lgwdp] + _b[L4.std_lgwdp] + _b[L5.std_lgwdp] + _b[L6.std_lgwdp] + _b[L7.std_lgwdp] + _b[L8.std_lgwdp]) / (1 - (_b[L1.std_mobility]))


/** unconstrained IDL model
use "Yourpath\080224_combineddata_clean.dta" 
save "Yourpath\092624_constrainedIDL.dta", replace

* Load the dataset
use "Yourpath\092624_constrainedIDL.dta", clear


*Get within R-square
* Identify unique state codes
levelsof fips_code, local(states)

* Loop through each state and run the unconstrained model
foreach state in `states' {
    
    di "Running model for state: `state'"

    * Preserve the full dataset before filtering
    preserve
    
    * Restrict data to the current state
    keep if fips_code == `state'
    
    * Check if there are any observations for this state
    count
    if r(N) == 0 {
        di "No observations for state: `state'. Skipping..."
        restore
        continue
    }

    * Standardize mobility and death variables
    egen zmobility = std(mobility)
    egen zlgwdp = std(lgwdp)

    * Generate mean-adjusted variables
    by fips_code: egen mzmobility = mean(zmobility)
    by fips_code: egen mzdeath = mean(zlgwdp)

    * Create detrended versions of mobility and death
    gen dzmobility = zmobility - mzmobility
    gen dzdeath = zlgwdp - mzdeath
set trace on
    * Generate lag variables for deaths (zlgwdp) for each state
    foreach i of numlist 1/20 {
        by fips_code: gen lag`i'_zlgwdp = zlgwdp[_n-`i']
        by fips_code: egen lag`i'_mzlgwdp = mean(lag`i'_zlgwdp)
        gen lag`i'_dzdeath = lag`i'_zlgwdp - lag`i'_mzlgwdp
    }
set trace off
    * Generate lagged mobility variables for each state
    foreach i of numlist 1/4 {
        by fips_code: gen lag`i'_zmobility = zmobility[_n-`i']
    }

    * Compute mean-adjusted lagged mobility
    by fips_code: egen mzmobility2 = mean(zmobility)
    gen dzmobility2 = zmobility - mzmobility2
	
* Check for missing values in key variables before regression
misstable summarize zmobility lag1_zlgwdp-lag20_zlgwdp

    * Estimate the Unconstrained IDL model using fixed effects for each state
    xtreg zmobility lag1_zlgwdp-lag20_zlgwdp spring fall summer i.cyclic, fe
    
    * Save the results to a file for the current state
    estimates store unconstrained_state_`state'
    esttab unconstrained_state_`state' using "results_state_`state'.txt", replace

    * Restore the dataset for the next state
    restore
}

* Combine all results into one log file for review
log using unconstrained_IDL_results.log, replace
foreach state in `states' {
    appendfile using "results_state_`state'.txt"
}
log close

set trace off


 *Export results to a .txt file
esttab results_*, mtitle unstack using "state_regression_results.txt"

 *Export results to a .csv file
esttab results_*, mtitle unstack using "state_regression_results.csv", csv
*/